Correspondence between Andreev reflection and Klein tunneling in bipolar graphene 
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Andreev reflection at a superconductor and Klein tunneling through an n-p junction in graphene 
are two processes that couple electrons to holes — the former through the superconducting pair 
potential A and the latter through the electrostatic potential U . We derive that the energy spectra 
in the two systems are identical, at low energies e ^ A and for an antisymmetric potential profile 
U{—x, y) = —U{x,y). This correspondence implies that bipolar junctions in graphene may have zero 
density of states at the Fermi level and carry a current in equilibrium, analogously to superconduct- 
ing Josephson junctions. It also implies that nonelectronic systems with the same band structure 
as graphene, such as honeycomb-lattice photonic crystals, can exhibit pseudo-superconducting be- 
havior. 

PACS numbers: 73.23.Ad, 73.23.Ra, 73.40.Lq, 74.45. +c 



Tunneling through an ri^ junction in graphene is 
called Klein tunneling [l|, d, y] with reference to relativis- 
tic quantum mechanics, where it represents the tunneling 
of a particle into the Dirac sea of antiparticles [4,] ■ Klein 
tunneling in graphene (see Fig. [1]) is the tunneling of an 
electron from the conduction band into hole states from 
the valence band — which plays the role of the Dirac sea. 
Several recent experiments ifl, 01 have investigated this 
unusual coupling of electron-like and hole-like dynamics. 

In the course of an analysis of these experiments a curi- 
ous similarity was noticed Q between negative refraction 
d, Q at an n-p junction and Andreev retroreflection 
at the interface between a normal metal (N) and a super- 
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FIG. 1: Conical band structure in graphene at two sides of a 
potential step (height 2/i, width d), forming an n-p junction. 
In equilibrium, all states below the Fermi level (indicated in 
blue) are filled and all states above are empty. Klein tunneling 
is the interband tunneling of an electron from the conduction 
band in the n region (blue ball at the right) into the valence 
band of the p region (blue ball at the left). In this work we 
show that the low-energy excitation spectrum of a symmetric 
n-p junction is the same as that of an NS junction, obtained 
by replacing the region x < by a superconductor. 
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FIG. 2: Periodic orbits in an n-p junction (panel a) and in a 
normal-superconductor (NS) junction (panel b), at e = in 
the case of an abrupt interface. (Solid and dashed lines dis- 
tinguish electron-like and hole-like trajectories.) Negative re- 
fraction in the n-p junction maps onto Andreev retrorefiection 
in the NS junction upon mirroring in the interface at a; = 0. 
Destructive interference of the electron-like and hole-like seg- 
ments of the periodic orbit suppresses the density of states at 
the Fermi level. Panels c and d show alternative geometries 
that exhibit a suppression of the local density of states in an 
unbounded system. 



conductor (S). As illustrated in Fig.[2^,b, the trajectories 
at an n-p junction and at an NS junction are related by 
mirroring x ^ —x at the interface (taken at x = 0). Here 
we show that the similarity is not limited to classical tra- 
jectories, but extends to the fully quantum mechanical 
wave functions and energy spectra. This implies that 
quantum effects associated with superconductivity, such 
as the proximity effect and the Josephson effect, have 
analogues in an n-p junction. 

We have found a precise mapping between the Dirac 
Hamiltonian (llj of an n-p junction and the Dirac- 
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Bogoliubov-De Gennes Hamiltoiiian [l2| of an NS junc- 
tion under the condition that the electrostatic poten- 
tial U in the n-p junction is antisymmetric, U{—x, y) = 
—U{x,y), with respect to the interface. The Fermi level 
is chosen at zero energy, symmetrically between the n 
and p regions. Such a symmetric n-p junction turns out 
to have the same excitation spectrum as an NS junction 
for excitation energies e small compared to the super- 
conducting gap A. After presenting the mapping in its 
mathematical form, we consider the two major physical 
implications: Zero density of states at the Fermi level 
and persistent current flow in equilibrium. A compari- 
son with computer simulations of a tight-binding model 
of graphene is presented at the end of the paper. 

Derivation of the mapping. — The correspondence be- 
tween Klein tunneling and Andreev reflection consists of 
a mapping of an eigenstate 5* of the Dirac Hamiltonian 
iJ of a symmetric n-p junction onto electron and hole 
eigenstates and '^h in the normal part a; > of the 
NS junction. The Dirac Hamiltonian is given (in the 
valley-isotropic representation) by 

H ^v[{p + eA) ■ (t](E}tq + U(To(E)To, (1) 

with p = —ih(d/dx,d/dy) the momentum operator in 
the X — y plane of the graphene layer, A = Bxy the vec- 
tor potential of a perpendicular magnetic field B, and v 
the electron velocity. The Pauli matrices Ui and act, 
respectively, on the sublattice and valley degree of free- 
dom (with CTo and tq a 2 x 2 unit matrix) . We introduce 
the time- reversal operator T = — (ctj, Ty)C, with C the 
operator of complex conjugation, and the parity opera- 
tor P = i{a-x (g) To)Tl, with TZ the operator of reflection 
(x I— > —x). The key property of the Dirac Hamiltonian 
that we need, in order to map the symmetric n-p junction 
onto an NS junction, is the anticommutation relation 

TVH = -HTP, (2) 

satisfied for any B when U{—x, y) ~ —U{x, y). 

Starting from a solution H"^ = e^* of the Dirac equa- 
tion in the n-p junction we now construct an eigenstate 
in the NS junction at the same eigenvalue e by means of 
the transformation 

■^e{x,y) = -^{x,y), ^!h{x,y)^V^{x,y). (3) 

According to Refs. [12, [3 the electron and hole wave 
functions ^'e, '^h iu the normal part of the NS junction 
should satisfy 

iJ^'e-e'J'e, -THT-^h^e-^h. 2; > 0, (4) 

with a boundary condition at the NS interface that for 
jej ^ A takes the form 

^-,,(0,^) = z((7, ® To)^'e(0,2/) = V^e{0.y)- (5) 

The proof of the mapping now follows by inspection: 
Firstly, Eq. ^ results directly from the transformation 



([3]) with the anticommutation relation ([2]). Secondly, 
since is continuous at a: = 0, the boundary condition 
([5]) is automatically satisfied. 

The applicability of the mapping extends to the crys- 
tallographic edges of the graphene layer in the following 
way: The edges of the n-p junction are described by the 
boundary condition "^{r) — M{r)^>{r) for r at the edge 
14, 15]. The mapping to an NS junction still holds, pro- 
vided that M commutes with V ^ which requires 

{a^®To)M{x,y) = M{-x,y)(a.^ (g, tq). (6) 

For example, an armchair edge parallel to the x-axis 
(with M (X Ux independent of x) satisfies the requirement 
([6]), but a zigzag edge parallel to the x-axis (M cx (Jz) 
does not. A pair of zigzag edges at x = zLW [with 
M{zLW,y) = ±az ® Tz], on the other hand, do satisfy 
the requirement An infinite mass boundary condi- 
tion [with M{±W, y) = ±(Ty ® t^], likewise, satisfies this 
requirement. 

Suppression of the density of states. — We have cal- 
culated the density of states p{e) by solving the Dirac 
equation in the n-p junction of Fig. [2^. The Fermi level 
(taken at e = 0) is separated from the Dirac point by the 
energy ±fj, in the n and p regions. We take an abrupt in- 
terface (width d small compared to the Fermi wave length 
Xp = hv/ii) and wide and long n and p regions (width 
W ^ Xp, length L 3> W). The precise choice of bound- 
ary condition at x = ±W does not matter in this regime, 
as long as it preserves the symmetry of the geometry. 

The calculation for the bipolar junction follows step- 
by-step the analogous calculation for the Josephson junc- 
tion in Ref. . The dispersion relation (smoothed over 
rapid oscillations) is given by 

e™(g) = TrET{m + - {hvq/n)^, \e\ <s: fi, (7) 

with m — 0,±1,±2,... the mode index and hq the 
momentum parallel to the n-p interface. (The energy 
Et — h.v/2W is the Thouless energy, which is ^ /i 
for W ^ Xp-) The resulting density of states p{e) = 
(4/7r) |9em/9(7|~^ is plotted in Fig. [31 It vanishes 
linearly as 

p(e)^Po\e\/ET (8) 

for small |e|, with po = (2p/Tr){hv)^'^ the density of states 
(per unit area and including spin plus valley degenera- 
cies) in the separate n and p regions. This suppression of 
the density of states at the Fermi level by a factor e / Et is 
precisely analogous to an NS junction, where the density 
of states is suppressed by the superconducting proximity 
effect (compare, for example, our Fig.[3]with Fig. 8 of Ref. 
[l^). In particular, the peaks in p{e) at £ = TrET{m-\- i) 
are analogous to the De Gennes-Saint James resonances 
in Josephson junctions [T7| . 

In a semiclassical description, the suppression of the 
density of states in the n-p junction can be understood 
as destructive interference of the electron-like and hole- 
like segments of a periodic orbit (solid and dashed lines in 
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FIG. 3: Density of states in the n-p junction of Fig. [2^, calcu- 
lated from Eq. ([7]). The dotted line is the value in the isolated 
n and p regions, which is energy independent for |e| <g /i. The 
density of states vanishes at the Fermi level (e = 0), according 
to Eq. dH). 
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FIG. 4: Persistent current through a ring containing an 
abrupt n-p interface, as a function of the magnetic flux 
through the ring. The solid curve is for zero temperature 
T = Q, the dashed curve for T — Er/'^kB, and the dotted 
curve for T = Er/ikB. 



Fig. [2^) . At the Fermi level, the dynamical phase shift 
accumulated in the n and p regions cancels, and what 
remains is a Berry phase shift of tt from the rotation of 
the pseudospin of a Dirac fermion [ll, [l^ . 

Persistent current. — If the n and p regions enclose a 
magnetic flux $, as in the ring geometry of Fig. [5] (inset), 
then the Berry phase shift can be compensated and the 
suppression of the density of states can be eliminated. 
The resulting flux dependence of the ground state energy 
E ~ A p{e)ede (with A the joint area of the n and 
p regions) implies that a current / — dE/d^ will flow 
through the ring in equilibrium at zero temperature, as 
in a Josephson junction According to Eq. ([H]), the 

order of magnitude 

lo = {e/h)El/5 = {e/h)NET (9) 

of this persistent current is set by the level spacing 6 = 
{Apo)^^ and by the Thouless energy Et — fiv/nr = 
Nd in the ring geometry (of radius r and width w <C 
r, supporting N = Ap.w/'Khv ^ 1 propagating modes). 
Because of the macroscopic suppression of the density of 
states, this is a macroscopic current — larger by a factor 
N than the mesoscopic persistent current in a ballistic 
metal ring 

We have calculated /(<&) for a simple model of an 
abrupt n-p junction in an iV-mode ring without inter- 
mode scattering, neglecting the effect of the curvature of 
the ring on the spectrum and also assuming that the 
magnetic field is confined to the interior of the ring. 
(These approximations are reasonable for Xp w r.) 
The slowly converging, oscillatory integral over p{e) was 
converted into a rapidly decaying sum over Matsub- 
ara frequencies by the method of Ref. [l^. The zero- 
temperature result is plotted in Fig. |4] (solid curve). The 
maximal persistent current is Ic ~ 0.2 /q. This is the 
same value, up to a numerical coefficient, as the critical 
current of a ballistic Josephson junction [2^. We have 



also included the results at finite temperature, showing 
the decay when the thermal energy ksT c± Et- 

Comparison with computer simulations. — To test 
our analytical predictions, we have performed computer 
simulations of a tight-binding Hamiltonian on a honey- 
comb lattice (lattice constant a). We took a symmetric 
n-p junction with zigzag boundaries at a; = ±W^ (with 
W/a = 400) and calculated the density of states p(e), 
smoothed by a Lorentzian (width 0.01 Et) to eliminate 
the rapid oscillations. Results are shown in Fig.O for dif- 
ferent Fermi wave lengths Xp = hv/ p, and widths d of the 
n-p interface [potential profile U{x) = — /itanh(4x/c?)]. 
A clear suppression of p{e) is observed within an energy 
range Et from the Fermi level at e = 0. The suppres- 
sion is somewhat smaller than predicted by Eq. ([5]) (black 
solid line), in particular for d ~ a (red curve, when the 
Dirac equation no longer applies) and for d > Af (blue 
curve, when Klein tunneling happens only near normal 
incidence Q). 

As expected, the suppression is sensitive to perturba- 
tions of the reflection symmetry. For example, as shown 
in Fig. [5] (yellow curve), a displacement of the n-p inter- 
face by Ai? spoils the systematic destructive interference 
due to the Berry phase, and thus eliminates the suppres- 
sion of the global density of states. 

We would still expect an effect on the local density of 
states if we could confine the carriers to the n-p interface. 
This might be achieved by means of the saddle point po- 
tential U = psign{xy) of Fig. [^t, or by means of the 
nonuniform magnetic field B = Bqx of Fig. [^Ji. Destruc- 
tive interference of the periodic orbits in each of these 
unbounded geometries will suppress the local density of 
states near the interface by the same mechanism as in 
the confined geometry of Fig. Because of disorder, 
the suppression will be limited to a mean free path or 
corrugation length from the n-p interface. Since the pre- 
dicted suppression of the density of states at the Fermi 
level happens at a large energy separation p from the 
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FIG. 5: Same as Fig. [3l but now calculated from a tight- 
binding model of graphene (lattice constant a, W/a = 400). 
The colors distinguish different values of Af and d, corre- 
sponding to an abrupt interface {Xp/a — 65, d/a — 12), a 
smooth interface {Xp/a = 12, d/a = 12), and an atomically 
sharp interface {Xp/a = 12, d/a ~ 1). The suppression of 
the density of states vanishes if the reflection symmetry is 
broken by displacing the interface (yellow curve, Xp/a = 65, 
d/a = 12, displacement — 65a). 



Dirac point (see Fig. [T]) , it should be distinguishable in a 
local measurement (for example, by a tunneling probe) 
from any features associated with the conical singularity 
in the band structure at the Dirac point. 

From a different perspective, the correspondence de- 
rived here offers the intriguing opportunity to observe 
superconducting analogies in non-electronic systems gov- 
erned by the same Dirac equation as graphene. An ex- 
ample would be a two-dimensional photonic crystal on 
a honeycomb or triangular lattice [3i llBl , in which the 
analogue of an n-p junction has been proposed recently 
[26l |. The detrimental effects of disorder should be rela- 
tively easy to avoid in such a metamaterial. 

We acknowledge discussions with J. Nilsson and R. A. 
Sepkhanov. This research was supported by the Dutch 
Science Foundation NWO/FOM. 



APPENDIX A: CALCULATION OF THE 
PERSISTENT CURRENT AND COMPARISON 
WITH SUPERCURRENT 

In this Appendix we present the calculation leading to 
the persistent current through the bipolar junction plot- 
ted in Fig.m We follow closely the analogous calculation 
for the supercurrent through a Josephson junction of Ref. 
[2^ . and compare the two systems at the end. For the 
sake of this comparison, it is convenient to work with 
the density of states p = {A/2)p per spin direction, inte- 
grated over the area A of the system. We will likewise, in 
this Appendix, count the number of propagating modes 
N = N/2 per spin direction. 



1. Persistent current 

The persistent current / — dF/d^ at temperature T is 
given by the derivative of the free energy F with respect 
to the flux $ enclosed by the ring containing the n-p 
junction. This can be expressed as an integral over the 
density of states, 

d r°° 

I^~2kBT— dep(e)ln[2cosh(£/2fcsT)]. (Al) 

'-'''^ J —ao 

We have set the Fermi energy at zero and used the 
electron-hole symmetry /5(e) — p{—s). The factor of two 
in front accounts for the two spin directions (which are 
not counted separately in p). 

Since the spectrum of the ring is discrete, the density 
of states p{e) = 5{e — Si) consists of delta functions 
at the solutions of the equation 

Tie) = To{£)Y[{£ - £'0 = 0. (A2) 

i 

(The index i counts spin-degenerate levels once.) The 
function ^0 is > and even in e, but can otherwise be 
freely chosen. The density of states is then written as 

p{e) = ---^Im\nT{e + iO+), (A3) 
TT de 

with 0"*" a positive infinitesimal. 

Substitution of Eq. into Eq. gives, using 

again the electron-hole symmetry, 

9k T d fOo+iO^ J 

/ ^ — / de ln[2cosh(£/2fcBT)]4-lnJ^(e). 

m d<P J„oo+io+ «e 

(A4) 

The expression for the persistent current becomes, upon 
partial integration, 

/ = — / tanh(e/2fcBT)lnJ="(e). (A5) 

We close the contour in the upper half of the complex 
plane. We assume that JTq is chosen such that In !F has no 
singularities for Ime > 0. The only poles of the integrand 
in Eq. (|A5p then come from the hyperbolic tangent, at the 
Matsubara frequencies itOn = (2n + l)iTTkBT. Summing 
over the residues we arrive at the expression [2^ 

d °° 

/ = -4fcBT— J2 In^(i^n)- (A6) 

n=0 

In our model of an A^-mode ring without intermode 
scattering we can calculate separately the contribution 
to / from each propagating mode, with transverse mo- 
mentum Qm- The total current is then a sum over these 
contributions, 

JV , 00 

/=-4fcBT^ — ^ln^(zu;„,g„). (A7) 

m— 1 n— 
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The function J^(e,q), which determines the energy levels 
in the bipolar junction for a given transverse mode, is the 
limit A ^ oo of the analogous function in a Josephson 
jimction fll|. We find 

J^ie,q) = ^2 sm6'+sm0_ 

+ cos6'+cos6'_ +cos(e$/?i), (A8) 
9± = Ej^^y/{n±e)'^ - {hvqy. (A9) 

Substitution into Eq. (jA7p gives the persistent current, 

N oo 

h 



I = ^keT^ sin(e$/?i) ^ ^ ^ 



t—l n—Q ^ ' J / 



(AlO) 



We wish to evaluate the expression (jAlOp in the regime 
fi :S> Et, N ^ 1. The sum over modes may be replaced 

by an integral, according to X]m=i ~^ i^/kp) Jq'' dq, 
with kp = fj,/hv the Fermi wave vector. Since the sum 
over the Matsubara frequencies converges exponentially 
fast for ujn ^ Et, we can also assume /i ;2> Ci;„. In this 
large- regime we may approximate 6± a ± ifi„, with 

a = Wi?T)[i-(g/MY/', (All) 

{ton/ Et)[1- (q/kpf ]-'^\ (A12) 

The function J- takes the form 

T{iujri,q) = X cos^ a + Y sin^ a, (A13) 
X ^ Z sinh^ rin + cosh^ fl^ + cos(e<I> /h), ( A14) 

Y = Zcosh^ il„ + sinh^ fi„ + cos(e$/?i), (A15) 



1 + {q/ki. 



^2 _ {Tivqf + {Ernnf 1 - {q/kpf ■ 



(A16) 



The phase a varies rapidly as a function of q, so we 
average 1/T first over this phase, 



1 



1 



da 

27r X cos^ a + Y sin'^ a 



1 

XY' 



(A17) 



We substitute Eq. (|A17j) into Eq. (jAlOp and evaluate it 
numerically, to arrive at the curves of / versus $ shown 
in Fig. m 



2. Comparison with supercurrent 

The mapping between bipolar junctions and Josephson 
junctions is illustrated in Fig.[6l Instead of a ring geome- 
try we may equivalently consider a planar SNS junction, 
with a phase difference (j) between the two superconduct- 
ing reservoirs. In the absence of mode mixing the two 
geometries carry the same supercurrent Ij, at the same 
number of transverse modes N. (The Thouless energy in 
the SNS junction is Et = Tiv/L, with L the separation 
of the two NS interfaces.) 




;Aoe-*/'^::;:;:;;;;;:Ao^ 




FIG. 6: Mapping of a bipolar ring containing two n-p junc- 
tions (left panel) onto a Josephson ring containing two NS 
junctions (right panels), by mirroring the hole- like trajecto- 
ries (dashed) in the line through the interfaces. The per- 
sistent current / through the bipolar ring at the left maps 
onto a supercurrent Ij through a Josephson ring at the right. 
Because the enclosed flux <I> is halved by the mapping, the 
ft/e-periodicity of I maps onto an h/2e periodicity of Ij. The 
flux enclosed by the Josephson ring in the upper right panel 
may be gauged away, with the introduction of a phase differ- 
ence (j> ~ e^/Ti between the order parameters at the two NS 
interfaces (lower right panel). 



The supercurrent Ij through the Josephson junction 
is given by [23| 



depj(e) ln[2 cosh(e/2fcBr)], 



(A18) 

with (j) the phase difference across the junction and pj 
the density of states per spin direction. The mapping 
relates (f> ^ e^/h (see Fig. (5]) and pj ^ p. Comparison 
of Eqs. (|XT|) and ((XT8|l then shows that ^ /($). 

The bipolar junction and Josephson junction therefore 
carry the same current in equilibrium. 

The result in the literature [H, [2^, [13, HH for a balhstic 
SNS junction is a piecewise linear dependence of Ij on 
(j) at zero temperature, close to but not identical to the 
solid curve in Fig.[4l As we will now show, the difference 
is due to the presence or absence of a step in the Fermi 
energy at the NS interfaces. 

On the one hand, the mapping between bipolar and 
Josephson junctions relies on the boundary condition ([5]) 
at the NS interface, which assumes that the Fermi energy 
fis in the superconductor is much larger than the value 
II in the normal region [T3|. On the other hand, Refs. 
[H, |2^ [13, HH assume /is = p,. The function T{e, q) is 
then given by 



T{e, q) — cos( 



(A19) 
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-n -n/l n/2 n 

FIG. 7: Supercurrent through a balhstic Josephson junction 
as a function of the phase difference (ji between the two super- 
conducting reservoirs, calculated from Eq. (|A20p . The solid 
curve is for zero temperature T = 0, the dashed curve for 
T = Er/AkB, and the dotted curve for T = Er/ikB. The 
difference with the analogous result for a bipolar junction in 
Fig. m arises because this figure is for equal Fermi energy 
/is = in superconductor and normal metal, while Fig. |4] 
maps onto a Josephson junction with /is 2> /i. 



(The (/^-dependence is repeated periodically outside of the 
interval — tt < (j) < tt.) We thus recover the piecewise 
linear 0-dependence of the supercurrent [H, [23, [3l| . 

For N ^ \ the sum over modes may also be evaluated 
as an integral, X]m=i ^ i^/kp) J^^'' dq, with the result 



IM = ^—^, H<^. (A22) 



The critical current If. — ireN Et / ^ti = (7r/8)/o is about 
two times larger than the maximal persistent current 
Ic ~ 0.2 /q found in the bipolar junction, because of the 
absence of a step in the Fermi energy at the NS interfaces. 

Eq. (|A22p holds in a two-dimensional geometry. In 

three dimensions the sum over modes becomes X]m=i ~^ 
{2N/kp) /p " qdq, resulting in 



resulting in 



N 



sm ( 



cosh2f2„((7m) -f cost 



(A20) 



The resulting supercurrent is plotted in Fig. [71 for differ- 
ent temperatures. 

At T = the sum over n reduces to an integral, 
J2^=o ~^ {2TrkBT)~^ duj, which evaluates to 



2eET 



N 



Y,[l-{q^/kpf]' 



/2 



oo 



ni—1 

duj ■ 



sm ( 



cosh 2uj + cos < 

N 



^Y.[l-[q,JkpfY/\ \^\<^. (A21) 



m— 1 



/jW^'/'^;^, 101 <7r, (A23) 



in agreement with Refs. [2^, [SO]. (The numerical coeffi- 
cient in Ref. [HI is different.) In the one-dimensional case 
iV = 1 of a single spin-degenerate mode (group velocity 



L^group 



;[1 - {qi/kpYY^'^) we find instead 



/7(0) = 



101 < 



(A24) 



in agreement with Ref . ,31j (up to a factor of two, pre- 
sumably because Ref. [3l| does not account for the spin 
degeneracy of the mode). 
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